1. Field of the Invention
This invention relates generally to the field of geophysical prospecting. More particularly, the invention relates to the field of seismic data processing. Specifically, the invention is a method for seismic migration using explicit depth extrapolation operators with dynamically variable operator length.
2. Description of the Related Art
The use of three-dimensional (3D) seismic methods has resulted in increased drilling success in the oil and gas industry. However, 3D seismic methods are still computationally expensive. A crucial point in processing 3D seismic data is the migration step, both because of its 3D nature and the computational cost involved. The accuracy, stability, and efficiency of 3D migration are determined by the wavefield extrapolation technique employed. Thus, it is desirable that the algorithms employed in 3D migration be accurate, stable, and efficient themselves, especially when 3D pre-stack migration is considered.
A number of the algorithms used in two-dimensional (2D) depth migration have not been successful in the 3D case. For example, implicit finite-difference extrapolation methods have the advantage of being unconditionally stable, but the disadvantage of being difficult to extend to the 3D case. The most common 3D implicit method is based on the operator splitting into alternating direction components. The splitting errors that these methods exhibit in the 3D case are translated into non-circularly symmetric impulse responses, which become unacceptable for dips higher than 45°. The splitting methods have errors that depend significantly on reflector dip and azimuth and thus have problems with reflector positioning errors. Similarly, two-pass methods have problems handling lateral velocity variations.
In contrast, explicit extrapolation methods that approximate the extrapolation operator as a finite-length spatial filter are easily extended to the 3D case. The difficulty with explicit extrapolation is that the stability condition is not automatically fulfilled. The stability condition is that no amplitude at any frequency will grow exponentially with depth. Stability must be guaranteed by careful design of the extrapolation operators.
Three-dimensional seismic wavefields may be extrapolated in depth, one frequency at a time, by two-dimensional convolution with a circularly symmetric, frequency- and velocity-dependent operator. This depth extrapolation, performed for each frequency independently, lies at the heart of 3D finite difference depth migration. The computational efficiency of 3D depth migration depends directly on the efficiency of this depth extrapolation. For these techniques to yield reliable and interpretable results, the underlying wavefield extrapolation must propagate the waves through inhomogeneous media with a minimum of numerically induced distortion over a range of frequencies and angles of propagation.
Holberg, O., “Towards optimum one-way wave propagation”, Geophysical Prospecting, Vol. 36, 1988, pp. 99–114, first disclosed explicit depth extrapolation with optimized operators. Holberg proposes a technique, solely for 2D depth migration, by generalizing conventional finite-difference expressions in the frequency-space domain. This technique yields optimized spatial symmetric convolutional operators, whose coefficients can be pre-computed before migration and made accessible in tables. The ratio between the temporal frequency and the local velocity is used to determine the correct operator at each grid point during the downward continuation, by fitting their spatial frequency response to the desired phase shift response over a range of frequencies and angles of propagation to control numerical distortion. Holberg's technique can be made to handle lateral velocity variations. However, the method only applies to 2D migration.
Blacquiere, G., Debeye, H. W. J, Wapenaar, C. P. A., and Berkhout, A. J, “3-D table driven migration”, Geophysical Prospecting, Vol. 37, 1989, pp. 925–958, extended the method disclosed in Holberg (1988) to 3D migration. The wavefield extrapolation is performed in the space-frequency domain as a space-dependent spatial convolution with recursive Kirchhoff extrapolation operators based on phase shift operators. The optimized operators are pre-calculated and stored in a table for a range of wavenumbers. The extrapolation is performed recursively in the space domain, so both vertical and lateral velocity variations can be handled. The Blacquiere et al. method is accurate, but is a full 3D method and hence computationally expensive.
Hale, D., “3-D depth migration via McClellan transformations”, Geophysics, Vol. 56, No. 11 (November, 1991), pp. 1778–1785, introduced a more efficient 3D scheme based on the McClellan transform, which gives numerically isotropic extrapolation operators. This commonly referred to as the Hale-McClellan scheme. Given the coefficients of one-dimensional frequency- and velocity-dependent filters similar to those used to accomplish 2D depth migration, McClellan transformations lead to an algorithm for 3D depth migration. Because the coefficients of two-dimensional depth extrapolation filters are never explicitly computed or stored, only the coefficients of the corresponding one-dimensional filters are required. Applying the two-dimensional extrapolation filter increases only linearly with the number of coefficients N in the corresponding one-dimensional filter, whereas the cost of convolution with a two-dimensional filter is generally proportional to N2. However, Hale's method has numerical anisotropy.
Two Soubaras publications: Soubaras, R, “Explicit 3-D migration using equiripple polynomial expansion and Laplacian synthesis”, 62nd Ann., Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1992, pp. 905–908, and Soubaras, R., “Explicit 3-D migration using equiripple polynomial expansion and Laplacian synthesis”, Geophysics, Vol. 61, No. 5 (September–October, 1996), pp. 1386–1393, improved the Hale-McClellan scheme. The Soubaras (1992, 1996) method used an expansion in second-order differential operators instead of the McClellan transform. Soubaras (1992, 1996) method also uses the Remez algorithm to design the coefficients for both the extrapolation operator and for the differential operators. This method avoids the numerical anisotropy to a large degree and is comparable in computational cost to the Hale-McClellan scheme.
The Soubaras approach takes advantage of operator circular symmetry, as do the McClellan transformations, but avoids the computation of a 2D filter approximating the cosine of the wavenumbers. The approach defines a Laplacian operator, approximates the Laplacian by the sum of two 1D filters approximating the second derivatives, and approximates the exact extrapolation operator by a polynomial. The second derivative operators and the polynomial expansion are both calculated by the Remez exchange algorithm.
Sollid, A., and Arntsen, B., “Cost effective 3D one-pass depth migration”, Geophysical Prospecting Vol. 42, 1994, pp. 755–716, makes the Soubaras (1992, 1996) scheme more cost effective. Sollid and Arntsen (1994) used frequency adaptive optimized derivative operators. The expansion in second order differential operators is used, but the expansion coefficients and differential operators are designed using a least squares approach rather than the Remez algorithm. A set of variable length second order differential operators for each wavenumber has different spectra and lengths to ensure that the resulting wave extrapolation is both accurate and efficient.
Mittet, R., 2002, “Explicit 3D depth migration with a constrained operator”, in 72nd Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, pp. 1148–1151, discloses a constrained explicit operator method which is a modification of the full 3D scheme discussed in Blacquiere et al. (1989), above, to make the scheme more computationally efficient. The number of independent operator coefficients is constrained to reduce the number of computer floating point operations required, thus increasing computer efficiency. The innermost coefficients in the core area of the extrapolation operator are computed in a standard fashion. The remaining outermost coefficients in the operator, related to very high dip and evanescent wave propagation, change only as a function of radius and are constant within radial intervals.
Thus, a need exists for an explicit depth extrapolation method for 2D and 3D seismic migration that is accurate, stable, and efficient.